* Replication of Figure 5, panel a, using Golder marginal effects plot
* You will need to specify the directory in which the replication data are stored in the second line of code 

clear
cd "[Directory with replication data here]"

version 11.0
#delimit ;
log using "sim5a.log", replace;
set more off;

use "fordham_wp2025.dta", clear ;

gen cotton100lbs=cotton*5;

logit multilat black time interactb ag interacta lynchings interactl cotton100lbs textapp, cluster(icpsr);

preserve;

set seed 339487731;

drawnorm SN_b1-SN_b10, n(10000) means(e(b)) cov(e(V)) clear;

postutil clear;
postfile mypost prob_hat0 lo0 hi0 prob_hat1 lo1 hi1 diff_hat diff_lo diff_hi 
            using "sim5a.dta", replace;
            
noisily display "start";

local a=0 ;
while `a' <= 26 { ;

{;

scalar h_time=17;
scalar h_cottonlbs=1.3;
scalar h_textapp=4.7;
scalar h_black=20;
scalar h_ag=20;
scalar h_lynchings=30;
scalar h_constant=1;

    generate x_betahat0  = SN_b1*h_black 
                            + SN_b2*h_time  
                            + SN_b3*h_black*h_time
                            + SN_b4*h_ag 
                            + SN_b5*h_ag*h_time
                            + SN_b6*h_lynchings 
                            + SN_b7*h_lynchings*h_time
                            + SN_b8*`a'
                            + SN_b9*h_textapp
                            + SN_b10*h_constant;
                            
    generate x_betahat1  = SN_b1*h_black
                            + SN_b2*(h_time+6)
                            + SN_b3*h_black*(h_time+6) 
                            + SN_b4*h_ag 
                            + SN_b5*h_ag*(h_time+6) 
                            + SN_b6*h_lynchings 
                            + SN_b7*h_lynchings*(h_time+6) 
                            + SN_b8*`a' 
                            + SN_b9*h_textapp
                            + SN_b10*h_constant;
    
    gen prob0 =normal(x_betahat0);
    gen prob1=normal(x_betahat1);
    gen diff=prob1-prob0;
    
    egen probhat0 =mean(prob0);
    egen probhat1=mean(prob1);
    egen diffhat=mean(diff);
    
    tempname prob_hat0 lo0 hi0 prob_hat1 lo1 hi1 diff_hat diff_lo diff_hi;  

    _pctile prob0, p(2.5,97.5) ;
    scalar `lo0' = r(r1);
    scalar `hi0' = r(r2);  
    
    _pctile prob1, p(2.5,97.5);
    scalar `lo1'= r(r1);
    scalar `hi1'= r(r2);  
    
    _pctile diff, p(2.5,97.5);
    scalar `diff_lo'= r(r1);
    scalar `diff_hi'= r(r2);  
   
    scalar `prob_hat0'=probhat0;
    scalar `prob_hat1'=probhat1;
    scalar `diff_hat'=diffhat;
    
    post mypost (`prob_hat0') (`lo0') (`hi0') (`prob_hat1') (`lo1') (`hi1') 
                (`diff_hat') (`diff_lo') (`diff_hi');
   
    };      
        
    drop x_betahat0 x_betahat1 prob0 prob1 diff probhat0 probhat1 diffhat;
	
    local a=`a'+ 1;
    
    display "." _c;
    
} ;

display "";

postclose mypost;

restore;

merge using "sim5a.dta";

gen MV = _n-1;

replace  MV=. if _n>26;

gen MV2=MV/5;

graph twoway hist cotton, percent color(gs14) yaxis(2) bins(100)
        ||  line prob_hat0 MV2, clwidth(medium) clcolor(blue) clcolor(black)
        ||  line lo0 MV2, clpattern(dash) clwidth(thin) clcolor(black)
        ||  line hi0 MV2, clpattern(dash) clwidth(thin) clcolor(black)
        ||  ,   
            xlabel(0 1 2 3 4 5.13, nogrid labsize(3)) 
            ylabel(0 0.2 0.4 0.6 0.8 1, axis(1) nogrid labsize(3))
            ylabel(0 5 10 15 20 25 30 35 40 45 50, axis(2) nogrid labsize(3))
            yscale(noline alt) yscale(noline alt axis(2)) xscale(noline) legend(off)
            yline(0, lcolor(black))
            yline(0 0.2 0.4 0.6 0.8 1, lcolor(white))
            xtitle("Cotton Production in 500-pound Bales Per Capita", size(3))
            ytitle("Probability of Voting in Support of Multilateralism", size(3))
            ytitle("Percent of Observations" , axis(2) size(3))
            xsca(titlegap(4)) ysca(titlegap(4)) ysca(axis(2) titlegap(4))
            scheme(s2mono) graphregion(fcolor(white) ilcolor(white) lcolor(white));
            
log close;
